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ABSTRACT 

In this paper we show how a Lagrangian variational principle can be used to derive the 
SPMHD (smoothed particle magnetohydrodynamics) equations for ideal MHD. We also con- 
sider the effect of a variable smoothing length in the SPH kernels after which we demonstrate 
by numerical tests that the consistent treatment of terms relating to the gradient of the smooth- 
ing length in the SPMHD equations significantly improves the accuracy of the algorithm. Our 
results complement those obtained in a companion paper rtPrice & M onaghan l2003l paper I) 
for non ideal MHD where artificial dissipative terms were included to handle shocks. 
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1 INTRODUCTION 

An advantage of deriving numerical algorithms from a variational principle is that conservation laws can be guaranteed. Another advantage 
is that the algorithms derived from a variational principle are often more stable than other algorithms. For example, in the case of smoothed 
particle hydrodynamics (SPH, for a review see Mona ghanll992l) . the density may be determined from the continuity equation, and it proves 
important for stability to combine the SPH continuity equation with the variational principle to deduce equations of motion. We call such a 
procedure consistent. 

bonet & Lokl ll999h have derived consistent SPH equations for fluids even when non standard forms of the continuity equation are used. 
They include the continuity equation as a constraint on Lagrangian density variations. The resulting equations possess very good stability 
properties when two fluids with very different densities, for example air and water, are in contact. Other, non consistent, forms of the SPH 
algorithm, for example with a standard acceleration equation but non standard continuity equation, exhibit instabilities. 

In the present paper we show how a Lagrangian variational principle can be used to derive the SPMHD (smoothed particle magneto- 
hydrodynamics) equations for ideal MHD. Variational equations for continuum MHD have been derived by Newcomb 1 1962) for both the 
Lagrangian and the Eulerian form of the equations (see also iHenvevll 1982^ lOppeneetll 19841 and iFieldl 1 986l) . In the Lagrangian form of the 
equations Newcomb makes use of flux conservation to relate changes in the magnetic field to changes in surface elements. In the present 
case, where we consider SPH particles, it is not clear how to prescribe such surface elements in a unique way from the particle coordinates. 
Instead we make use of the induction equation in its Lagrangian form and treat this as a constraint. An alternative, which we do not explore 
here, is to begin with plasma physics and prescribe the fields in terms of currents. Such an approach would be natural for particle methods 
(e.g. PIC) which have been so effective for plasma physics where the electrons would be treated as one fluid and the ions as another. 

The plan of this paper is derive the equations of motion from a standard Lagrangian for SPH particles with either, or both, the continuity 
and induction equations treated as constraints (0. We then consider the effect of variable smoothing length in the SPH kernels (Sj4j after 
which we demonstrate by numerical tests that consistent treatment of the variable smoothing length in the SPH equations significantly 
impro ves the accuracy of SPMHD shocks and of wave propagation (J6j. Our results complement those obtained in a companion paper 
IPrice & Monagh arl2003l hereafter paper I) for non-ideal MHD where artificial dissipative terms were included to handle shocks. 



2 THE LAGRANGIAN 

Variational principles for MHD have been discussed by many authors (e. g. iNewcombl 1 962LlHenvevl 1 982LIOppeneeJ 1 984t IFieldl 1 9861) and 
the Lagrangian is given by 
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which is simply the kinetic minus the potential and magnetic energies. The SPH Lagrangian is therefore 

L S ph = y ' ) 



-Vfc - U b {p b ) - — — 

2 2 fio pb 



(1) 
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where we have replaced the integral with a summation and the volume eleme nt pdV wi th the mass per SPH particle m. Variational principles 
for SPH in relativistic and non-relativistic fluid dynamics have been given by M onaghan & Price! 1200 lh . 



3 SPMHD EQUATIONS 
3.1 Equations of motion 

3.1.1 Standard formulation 

Ideally we would wish to express all the terms in the Lagrangian in terms of the particle co-ordinates, which would automatically guarantee 
the co nservation of momentum and energy since the equations of motion result from the Euler-Lagrange equations (e.g. Monag han & Pricel 
l200ll) . The density can be written as a function of the particle coordinates using the usual SPH summation, that is 

(3) 



pa = 2_j mb Wab, 



where W a b = W(r a —Tt,h) is the SPH interpolation kernel. Taking the time derivative of this expression, we have the SPH version of the 
continuity equation 



dp a 

dt 



m b (v a - Vfc) ■ VaWab 



(4) 



The internal energy is regarded as a function of the density, where from the first law of thermodynamics we have 

du _ P 
dp p 2 

The magnetic field is evolved in SPH according to 

^ = — y^m b [B a (Vab-V a Wab)-Vab(~B a -V a Wab)], 

at p a ' 
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or equivalently 
B 

PJa * 

(e.g- IPhillips & Monaghanlll985tlMonaghaniri992l pa per I). We note that these equ ations represent the correct formulation of the induction 
equation even in the presence of magnetic monopoles (lanhunen 2000; Dellai 2001). 

However it is not intuitively obvious how the magnetic field B should be related to the particle co-ordinates, or even that it could be 
expressed in such a manner (in the SPH context this would imply an expression for B such that taking the time derivative gives ^6) or jTJ, 
analogous to for the density), though it could be done easily for a plasma with the electrons a nd ions described by separate sets of SPH 
particles. We may however proceed by introducing constraints on B in a manner similar to that of lBonet & Lokl Il999h . that is we require 



idt = / <5Ldt = 0, 



where we consider variations with respect to a small change in the particle co-ordinates 5r a . We therefore have 

2 
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The Lagrangian variations in density and magnetic field are given by 
Spb = m c (Sr b - 5r c ) ■ V& W bc 



(8) 



(9) 
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where we have used and Q respectively (note that we also recover the following results if we use instead of Q). Using {ToJ, {Tl) and 
in {5) and rearranging, we find 
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where S a b refers to the Kronecker delta. Putting this back into JSJ, integrating the velocity term by parts and simplifying (using V a W a b = 
—VbWba), we obtain 
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The SPH equations of motion are therefore given by 
dvl \ - 



dt 
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nib 
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where the stress tensor S 13 is defined as 



-PS tJ + 



Mo V 



(13) 



(14) 



(15) 



This form of the m agnetic force term conserves linear momentum exactly (angular momentum is discussed in but was shown 
bv IPhillips & Monaghanl (l985) to be unstable i n certain regimes ( low magnetic /3). We resolve this instability by adding a short range 
repulsive force to prevent particles from clumping lMonagharT (2000), the implementation of which is described in paper I. We note that the 
conservative form of the momentum equatio n was derived using a non-conservative induction equation, which agrees with the treatment of 
magnetic monopoles suggested bv ljanhunerJ kood) and lDellaJi200lll 



3.1.2 Alternative formulation 



Consistent sets of SPMHD equations may also be derived using alternative forms of the continuity and induction equations. We give the 
example below since alternative forms of the pressure terms in the momentum equation are often explored in the context of SPH, without 
alteration of the other equations to make the formalisms self-consistent. We expect that a lack of consistency in the discrete equations will 
inevitably lead to loss of accuracy in the resulting algorithm. For example, using the continuity equation 

dp a 



,. = Pa ^ m b -VaWab, 

dt t-^ p h 

b 



and the induction equation 



d_ 

dt 
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p^ 4^ 
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(16) 



(17) 



results in the momentum equation 
dv a _ sr^ 



dt 



b 



m b 



PaPb 



VlWab 



(18) 



This form of the SPMHD equations also conserves linear momentum exactly and in the hydrodynamic case has been found to give better 
performance in situations where there are large jumps in density (for example at a water-air interface). The consistent form of the energy 
equations are given in £13.2.31 
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3.2 Energy equation 

3.2.1 Internal energy 

The internal energy equation follows from the use of the first law of thermodynamics, that is 

dUq _ P a dp a 

dt p a dt 

Using the standard continuity equation J4j therefore gives 



du a P a \ -> ... 
—TT = —o > rn b V ab ■ VaWab 

dt pi ^ 



(19) 



(20) 



3.2.2 Total energy 

The Hamiltonian is given by 



(21) 



which represents the conserved total energy of the SPH particles since the Lagrangian does not explicitly depend on the time co-ordinate. 
Using we have 
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Taking the (comoving) time derivative, we have 
dE \ - 



dt 



v rfvq | dug dp a ^ 1 Bj dpa ^ ^ d /B a 
a dt dp a dt 2 p a dt a dt \ p a 



(22) 



(23) 



where the first term is specified by use of the momentum equation i 14t . the second term using the first law of thermodynamics |5J and the 
continuity equation J4j, the third term by the continuity equation J4j and the fourth term by the induction equation 0. Using these equations 
and simplifying we find 



dE 
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such that the total energy per particle is evolved according to 



de a 
dt 

where 



E 

b 



„2 



v b - 



1 2 
ea = -V a 



2 p a 



(24) 
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(26) 



is the energy per unit mass. 



3.2.3 Alternative formulation 

For the alternative formulation given in i|3.1.2l the internal energy equation is given by 

du a Pa V ab 



dt 



Pa 



E 



rrib- 



Pb 



VaWab, 



and the total energy equation by 



de a 
dt 



E 

b 



S>6 + S>q 
papb 



(27) 



(28) 



4 VARIABLE SMOOTHING LENGTH TERMS 



The smoothing length h determines the radius of interaction for each SPH particle. Early SPH simulations used a fixed smoothing length for 
all particles, however allowing each particl e to have its own assoc iated smoothing length which varies according to local conditions increases 
the spatial resolution substantially faernauist & KatJl989l:lB"enzll99d) . The usual rule is to take 
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/ 1 xdM 

h a oc ( — J , (29) 

where v is the number of spatial dimensions, although others are possible jMonaghaifcOQol Implementing this rule self-consistently is more 
complicated in SPH since the density p a i s itself a function of the smoothing length h a via the relation The usual rule is to take the time 
derivative of 1291 . giving (e.g. lBenzll990h 

dh a _ ha_dp_ 

dt vp a dt ' 

which can then be evolved alongside the other particle quantities. 

This rule works well for most practical purposes, and maintains the relation 1291 particularly well when the density is updated using 
the continuity equation J4j. However, it has been known for some time that, in orde r to be fully self-consistent, extra terms involving th e 
derivative of h should be included in the momentum and energ y equations (e.g.lNelso n|l994l: lNelson & Papaloizoill994l : ISerna et alll 9961. 
Attempts to do this were, however, complicated to implement l Nelson & Papaloizou 1994) and therefore not generally adopted by the SPH 
community. Recentlv |Springel & Hemquistl i2002t) h ave shown that the so-called V /i terms can be self-consistently included in the equations 
of motion and energy using a variational approach, jspringel & HernquisJ i2002t) included the variation of the smoothing length in their 
variational principle by use of Lagrange multipliers, however, in the context of the discussion given in j3] we note that b y expressing the 
smoothing length as a function of p we can therefore specify h as a function of the particle co-ordinates iMonaghan 2002). That is we have 
h = h(p) where p is given by 

pa = ^ m b W(r a b, h a )- (31) 

b 

Taking the time derivative, we obtain 

^ = ^-^m b w ab -V a W ab {h a ), (32) 

° b 



where 



a, 



■ . dh a ST^ dWab(h a ) 
Op a ^— ' dh a 



(33) 



The equations of motion in the hydrodynamic case may then be foun d using the Euler-Lagrange equations a nd will automatically 
conserve linear and angular momentum. The resulting equations are given by fcniwel & Hexnnuisl200i liv Tor,aph a rl2002h 



^ = -Vm 6 J\x7 a Wab(ha) + -j^j V aW a b(h b ) . (34) 
dt [Q a pi Sib Pi 

Calculation of the quantities Q, involve a summation over the particles and can be computed alongside the density summation 13 1 i . 
To be fully self-consistent <3 1 1 should be solved iteratively to determine both h and p self-consistently. In prac tice this means that an extra 
pass over the density summation only occurs when the density changes significantly between timesteps. lSpringel & Hernquistl 120021) also 
suggest using the continuity equation 1321 to obtain a better starting approximation for p and consequently h. We perform simple fixed point 
iterations of the density, using a predicted smoothing length from 1301 . Having calculated the density by summation, we then compute a new 
value of smoothing length h new using <29> . The convergence of each particle is then determined according to the criterion 

lhnew 7 - fe| < 1.0 x HT 2 . (35) 
h 

We then iterate until all particles are converged, although for efficiency we do not allow the scheme to continue iterating on the same 
particle(s). Note that a particle's smoothing length is only set equal to h new if the density is to be recalculated (this is to ensure that the 
same smoothing length that was used to calculate the density is used to compute the terms in the other SPMHD equations). The calculated 
gradient terms 1331 may also be used to implement an iteration scheme such as the Newton-Raphson method which converges faster than our 
simple fixed point iteration. We also note that in principle only the density on particles which have not converged need to be recomputed, 
since each particle's density is independent of the smoothing length of neighbouring particles. These considerations will be discussed further 
in the multidimensional context since the cost of iteration is of greater importance in this case. 

Since we cannot explicitly write the Lagrangian J5J as a function of the particle co-ordinates, we cannot explicitly derive the SPMHD 
equations incorporating a variable smoothing length. We may, however deduce the form of the terms which should be included by consistency 
arguments. We start with the SPH induction equation in the form 



(36) 



Expanding the left hand side, we have 
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— 7T = > m b v ab (B a ■ V a W ab ) (37) 

dt p a ^ pa dt 

b 

If the smoothing length is a given function of the density, then the SPH continuity equation is given by I32i and 1371 becomes 

-— Vm t (v a6 (B a ■ V a VM - ^-B a [v ai) • V„Wai(h.)]) ■ ( 38 ) 

Pa ' I "a J 



However in one dimension these terms must cancel to give B x — const, and thus we deduce that the correct form of the induction equation 
is therefore 



^ = -7^— V m b {Vab [Ba ■ VaWab(ha)] - B a [v a b ■ V aW a b(ha)]} . 



or in the form |36j we would have 

— 2 ^ m b V ab [B a • VaWab(/l a )]. 



d/B 



Flap a 



(39) 



(40) 



Using 1391 or 1401 and 1321 as constraints we may then derive the equations of motion using the variational principle described in j3| to gi ye 
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dt. ~ 
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The total energy equation is given by 



de a \ - 



Q.P 1 



77^ ) vlViWab(ha) +(77^) V a ViW ab {h 



n P ' 2 



5"' 
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(41) 
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We show in H6.1l that including the correction terms for a variable smoothing length in this manner significantly improves the numerical 
wave speed in the propagation of MHD waves and enables the shock tube problems considered in paper I to be computed with no smoothing 
of the initial conditions. 



whilst the internal energy equation is found using the first law of thermodynamics and 1321 . that is 

^7 = ^ m b Vab ■ VaWab(ha) 

dt Q, a pl ^ 



5 MOMENTUM CONSERVATION 

The equations of motion conserve linear momentum exactly. However, angular momentum is not conserved exactly because the stress force 
between a pair of particles is not along the line joining them. Returning to <41> . and considering motion in 2 dimensions x and y, the change 
in angular momentum of the system is given by 

-17 ^( r<l x V ") Z ~ y^y^ m a m b ([d-ab ~ all) VabX a b + ^tblVab - x 2 ab }) F ab , (44) 

a a b 

where y a b = ya — yb, x a b = x a — x b , u 1 -' = S 1 -' /Q,p 2 and a'l = a l J + a % b 3 . We have replaced VW a b by v a bF a b- It can be seen from J44I 
that if the stress is isotropic, and proportional to the identity tensor, as is the case for isotropic fluids, the rate of change of angular momentum 
vanishes. If, however, the stress is not proportional to the identity tensor then the total angular momentum of the system will change. It can 
be shown that when the summations can be converted to integrals the angular momentum is conserved exactly. 

The same problem arises in the case of elastic stresses where the problem is exacerbated by the fact that the particles near the edge of 
a solid have densities similar to the interio r and the particles do not have neighbours exterior to the solid. In this case the conservation of 
angular momentum is significantly in error, bonet & Lokl Il999h showed, however, that angular momentum could be conserved by altering 
the gradient of the kernel to a matrix operator. The astrophysical problem could be solved in the same way but we expect the astrophysical 
conservation to be very much better without changing the kernel, because edges are associated with low density and correspondingly low 
angular momentum. 



6 NUMERICAL TESTS 

We de monstrate the usefulness of the variable smoothing length terms in the MHD case by the simulation of MHD waves and the lBrio & Wul 
1 1988) shock tube problem. We find that with the variable smoothing length terms included it is better to use 1401 to update the magnetic 
field rather than 1391 since we find that using 1391 can lead to negative thermal energies in the shock tube problem. The results shown use the 
total energy equation (as in paper I) although the similar results are obtained when the thermal energy is integrated. 
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6.1 MHD waves 

The equations of magnetohydrodyna mics admit three 'families' of waves, the so called slow, alfven and fast waves (appendix 0. The tests 
presented here are taken from lDai & Woodward 1 199^) . We consider travelling slow and fast MHD waves propagating in a ID domain, where 
the velocity and magnetic field are allowed to vary in three dimensions. We use 7 = 5/3 for the problems considered here. The perturbation 
in density is applied by perturbing the particles from an initially uniform setup (since we use equal mass particles). Details of this perturbation 
are given in appendix |B]and the amplitudes for the other quantities are derived in appendix|A| We leave the artificial dissipation on for this 
problem, with the llVIorris & Monaghanl 1 19971) switch implemented using Kmin = 0.05 (see paper I for details of this implementation). This 
is to demonstrate that the switch is effective in turning off the artificial dissipation in the absence of shocks. The variable smoothing length 
terms (0 do not affect the wave profiles but inclusion of these terms gives very accurate numerical wave speeds. 

The fast wave is shown in FigureQ with the dashed line giving the initial conditions. The initial amplitude is 0.55% as in lPai & Woodwarc 

(1998). Results are shown at t=10 for five different simulations using 32, 64, 128, 256 and 512 particles in the x-direction. The properties 
of the gas are set such that the fast wave speed is very close to unity, meaning that the solution at t = 10 should be in phase with the initial 
conditions. FigureQdemonstrates that this is accurately reproduced by the SPMHD algorithm. The effects of the small amount of dissipation 
present can be seen in the amount of damping present i n the solutions. The sma ll amount of steepening observed wave profiles is due to 
nonlinear effects and agrees with the results presented bv lDai & Woodward (l998). 

The slow MHD wave is shown in Figure|5| again with the dashed line giving the initial conditions. The perturbation amplitude is 0.6% as 
in lDai & Woodward] Il998l) . Results are again shown att = 10 at resolutions of 32, 64, 128, 256 and 512 particles. In this case the properties 
of the gas being set such that the slow wave speed in the medium is very close to unity, again meaning that the solution at t = 10 should be 
in phase with the initial conditions. We see in Figure[2]that this is reproduced by the SPMHD solution for the higher resolution runs. The 
artificial dissipation is again turned on using the switch and a minimum of K min — 0.05. The wave is slightly overdamped in this case since 
we construct the artificial dissipation using the fastest wave speed (c.f. paper I) which in this case is approximately three times the wave 
propagation speed. 



6.2 Shock tube 

As an additional example o f the advantages of the consistent smoothing length evolution and the variable smoothing length terms we recal- 
culate the lBrio & wJ il988h shock tube test from paper I. In this case however we apply no smoothing whatsoever to the initial conditions 
and calculate the solution using the density summation <3H . the total energy equation 1421 and the induction equation I40i . As in paper I 
we set up the problem using approximately 800 equal mass particles in the domain x — [—0.5, 0.5]. Conditions to the left of the shock are 
given by (p, P, v x , v y , B y ) = [1, 1,0,0, 1] and to the right by (p, P,v x ,v y , B y ) = [0.125, 0.1, 0, 0, -1] with B x = 0.75 and 7 = 2.0. The 
results are shown in Figure|5|at time t = 0.1 and may be compared with the numerical solution from lBalsaraHl998T) given by the solid line. 
The results may also be compared with Figure 2 in paper I. The non-smoothed initial conditions result in a small starting error at the contact 
discontinuity and a small overshoot at the end of the rarefaction wave, however the compound wave in particular is significantly less spread 
out than in the results given in paper I. The consistent update of the smoothing length discussed in ^results in some extra iterations of the 
density (for most of the simulation two passes over the density summation are used). 

As a final example we also recompute the shock tube test shown in Figure 7 of paper I. The initial conditions to the left of the shock are 
given by (p, P, v x ,v y ,v z , B y , B z ) = [1,1, 36.87, —0.155, —0.0386, 4/ (47T") 1 / 2 , 1/ (47r) 1//2 ] and to the right by (p, P, v x ,v y ,v z , B y , B z ) = 
[1, 1, -36.87, 0,0, 4/(47r) 1/2 , l/(47r) 1/2 ] withB^ = 4.0/(47r) 1//2 and 7 = 5/3, resulting in two extremely strong fast shocks which propa- 
gate away from the origin. The resolution varies from 400 to 1286 particles throughout the simulation due to the inflow boundar y cond itions. 
Results are shown in Figure[4]at time t — 0.03 and compare extremely well with the exact solution given bv lDai & Woodward fl994h (solid 
lines). In paper I the post-shock density and transverse magnetic field components were observed to overshoot the exact solution. In Figure 
|4|we observe that these effects are no longer present when the variable smoothing length terms are self-consistently accounted for. 



7 SUMMARY 

In summary, we have shown that 

(i) The equations of motion and energy for SPMHP can be derived from a variational principle using the continuity and induction 
equations as constraints. This demonstrates that the equation set is consistent and the resulting equations conserve linear momentum and 
energy exactly. In the MHP case this also demonstrates that the treatment of source terms proportional to V ■ B is consistent, as discussed in 
paper I with reference to Janhunen 1 200C) and Pellar 1 2001). 

(ii) The correction terms for a variable smoothing length may be derived naturally from a variational approach. Accounting for these terms 
is shown to improve the accuracy of SPH wave propagation. 
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Figure 1. Results for the ID travelling fast wave problem. Initial conditions are indicated by the dashed line. Results are presented after 10 periods for 
simulations with 32, 64, 128, 256 and 512 particles. The fast wave speed in the gas is very close to unity which is accurately reproduced by the SPMHD 
solution (ie. the numerical solution is in phase with the initial conditions). The artificial dissipation is turned on but uses the switch of lMorris & Monaahan 
1 1997) which dramatically reduces its effects away from shocks. The wave is steepened slightly by nonlinear effects. 
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Figure 2. Results for the ID travelling slow wave problem. Initial conditions are indicated by the dashed line and results are presented after 10 periods for 
simulations with 32, 64, 128, 256 and 5 1 2 particles. The slow wave speed in the gas is very close to unity, such that the numerical solution at t = 10 should be 
in phase with the initial condit ions. This is we ll represented by the SPMHD solution for the higher resolution runs. The artificial dissipation is turned on but 
uses the switch of Morris & Monaehan 1 1997) which dramatically reduces its effects away from shocks. The wave is steepened slightly by nonlinear effects. 
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Figure 3. Results of the lBrio & Wul Il988i) shock tube test with no smoothing of the initial conditions. Initial conditions to the left of the origin are given 
by (p,P,v x ,v y ,By) = [1, 1, 0,0, 1] and to the right by (j>, P,v x ,v y ,B y ) = [0.125,0.1,0,0,-1] with B x = 0.75 and 7 = 2.0. Profiles of density 
pressu re, v x , v y , transverse magnetic field and thermal energy are shown at time t = 0.1 and may be compared with the numerical solution from Balsara 
1 1 998 ) given by the solid line. In this case the density summation, total energy equation and the induction equation using B/p have been used, incorporating 
the variable smoothing length terms. 
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Figure 4. Results of the MHD shock tube test with initial conditions to the left of the shock given by (p, P, v x , v y , v z , B y , B z ) = 
[1, 1,36.87, -0.155, -0.0386, 4/ (47T) 1 / 2 , l/(4?r) 1 / 2 ] and to the right by (p, P, v x ,v y ,v z , B y , B z ) = [1, 1, -36.87, 0, 0, 4/( 4tt) 1 / 2 , 1/(4tt) 1 / 2 1 with 
B x = 4.0/(47t) 1/ ' 2 and 7 = 5/3. Results are shown at time t = 0.03 compare extremely well with the exact solution given bv lDai & Woodwardl Il993) 
(solid lines). The overshoots in density, pressure and magnetic field observed in paper I are no longer present due to our self-consistent inclusion of terms 
relating to the gradient of the smoothing length. 
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APPENDIX A: LINEAR WAVES IN MHD 

In this section we describe the setup used for the MHD waves described in £16.11 The MHD equations in continuum form may be written as 

% = -pV-v, (Al) 

at 

dv VP B x (V x B) , A(Vl 

-77 = " -, (A2) 

at p nop 

^ = (B-V)v-B(V-v), (A3) 
together with the divergence constraint V ■ B = 0. We perturb according to 



p = 


po + 5p, 


V 


v, 


15 = 


Bo + SB, 


SP = 


c s 8p. 


where c 2 a 


= fPo 1 po is the sound speed. Considerinj 


d(Sp) 
dt 


= -po(V-v), 


dv 


2 V((Sp) B„ x (V x <5B) 


lit 


po p-oPo 


d(SB) 
dt 


= (Bo- V)v-Bo(V-v). 


Specifyin: 


g the perturbation according to 


Sp = 




v = 


ve *(kx—t) 5 


SB = 


be <(te-«*) ) 


we have 




-loD 


= -po(v-k) 


—UJV 


= -d— - — [(Bo-b)k-(Bo 

Po popo 


— coh 


= (Bo -k)v-Bo(k-v). 



(A4) 

(A5) 
(A6) 
(A7) 



(A8) 
(A9) 

>] (A10) 

(All) 

Considering only waves in the x-direction (ie. k = [k x , 0, 0]), defining the wave speed v = co/k and using <A9> to eliminate D, equation 
gives 

Byob y + B z0 b 2 



popo 



B x oby 
Popo 
B x ob z 



(A 12) 
(A 13) 

(A14) 

popo 

where b x — since V • B = 0. Using these in <A1 H we have 

Vby = — B X QVy + ByOVx, (A15) 

vb z = -BxoVz + B z0 v x - (A16) 

We can therefore solve for the perturbation amplitudes v x ,v y ,v z , b y and 6 Z in terms of the amplitude of the density perturbation D and the 
wave speed v. We find 

v x = — (A17) 
P 
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b z v 



B 2 \ 

^X \ 


J3x By 


pop) 


Pop 


bV\ 


B X B Z 


Pop J 


pop 


bV\ 
pop J 


— VByV X 


POp) 


= vB z v x 



where we have dropped the subscript 0. The wave speed v is found by eliminating these quantities from <A12> . giving 



(v 2 - B'i/fiop) 



v — V c» 



B x + B y + B z 
Pop 



+ ■ 



Pop 



0, 



(A18) 
(A 19) 
(A20) 
(A21) 

(A22) 



(A23) 



which reveals the three wave types in MHD. The Alfven waves are those with 

y 2 _ B 2 
Pop' 

These are transverse waves which travel along the field lines. The term in square brackets in <A22t gives a quartic for v (or a quadratic for 
v 2 ), with roots 



c s + 



B, + B,. + B, 



POP 



±1 ci + 



B 2 + B 2 +B 2 \ c 2 B 



Pop 



POP 



(A24) 



which are the fast(+) and slow(-) magnetosonic waves. 



APPENDIX B: DENSITY PERTURBATION IN SPH 

The perturbation in density is applied by perturbing the particles from an initially uniform setup. We consider the one dimensional perturba- 
tion 

p = po[l+Asin(ka:)], (Bl) 
where A = D/po is the perturbation amplitude. The cumulative total mass in the x direction is given by 

M(x) = p J [1 + Asin(kx)]dx 

= po[x- Acos(ka;)]S, (B2) 
such that the cumulative mass at any given point as a fraction of the total mass is given by 

Mix) {m) 

M(Xmax) ' 

For equal mass particles distributed in x = [0, x max ] the cumulative mass fraction at particle a is given by x a /x max such that the particle 
position may be calculated using 

_ M(Xq) 

- TT7Z V ( B4 ) 

Substituting the expression for M(x) we have the following equation for the particle position 

x a x a - Acos(ka: a ) _ Q 

Xmq X \Xmq X AcOs(k^ ma:c )] 

which we solve iteratively using a simple Newton-Raphson rootfinder. With the uniform particle distribution as the initial conditions this 
converges in one or two iterations. 
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